Relativistic orbits with gravitomagnetic corrections 
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Corrections to the relativistic orbits are studied considering higher order approximations induced 
by gravitomagnetic effects. We discuss in details how such corrections come out taking into account 
"magnetic" components in the weak field limit of gravitational field and then the theory of orbits 
is developed starting from the Newtonian one, the lowest order in the approximation. Finally, 
the orbital structure and the stability conditions are discussed giving numerical examples. Beside 
qq \ the standard periastron corrections of General Relativity, a new nutation effect is due to the c -3 

corrections. The transition to a chaotic behavior strictly depends on the initial conditions. The 
f^*) ■ orbital phase space portrait is discussed. 
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(N ; I. INTRODUCTION 

The analogy between the classical Newton and Coulomb laws led to investigate if masses in motion, considered as 
O" 1 charges, could give rise to a "gravitational" magnetic field. 

5^ ' In fact, the magnetic field is produced by the motion of electric-charge, i.e. the electric current. The analogy 
consists of the fact that a mass-energy current can produce what is called "gravitomagnetic" field. 

The pioneering approach to the problem is due to Maxwell himself which, in one of his fundamental works on 
electromagnetism, turned his attention on the possibility to formulate the theory of gravitation in a form corresponding 
J> ■ to the electromagnetic equations However, he was puzzled by the problem of the energy of the gravitational 
C*") \ field, i.e. the meaning and the origin of the negative energy due to the mutual attraction of material bodies. In 
fact, according to him, the energy of a given field had to be "essentially positive", but this is not the case for the 
gravitational field. To balance this negative energy, a great amount of positive energy is required, in the form of 
energy of the space (a sort of back-reaction). But, since he was unable to understand how this could be, he did not 
proceed further along this line of thinking since the problem can be addressed and solved only in the framework of 
General Relativity. 

Later, Holzmuller [2] and Tisserand [3| proposed to modify the Newton law introducing, in the radial component 
of the force, a term depending on the relative velocity of the two attracting particles (see also @, Hj]). Also Heaviside 
@, 0] investigated the analogy between gravitation and electromagnetism considering the propagation of gravitational 
energy in terms of a sort of gravito - electromagnetic Poynting vector: however, also in this case, he failed to frame 
the problem of gravitational energy in a self-consistent scheme. 

Finally, the formal analogy between electromagnetic and gravitational fields was explored by Einstein 0, in the 
framework of General Relativity, and then by Thirring [§]. This author pointed out that the geodesic equation can be 
written as a Lorentz force splitting the gravitational field in gravito-electric and gravito-magnetic components. The 
final result of these studies was that any theory which combines Newtonian gravity together with Lorentz invariance 
in a consistent way, has to include a gravitomagnetic field, which is generated by the mass-energy current. This is 
the case, of course, of General Relativity: it was shown by Lense and Thirring [l(|, that a rotating mass generates 
a gravitomagnetic field, which in turn, causes the precession of planetary orbits. To be more precise, H. Pfister has 
recently shown that it would be better to speak about an Einstein - Lense - Thirring effect 

It is interesting to notice that also Lodge and Larmor, at the end of the nineteenth century, discussed the effects of 
frame dragging on a non- rotating interferometer [12], but within the framework of an aether-theoretical model. This 
frame dragging corresponded, in fact, to the Lense-Thirring effect of General Relativity. However, at the beginning 
of the XX century, when Lense and Thirring published their papers, the effect named after them, which is indeed 
very small in the terrestrial environment, was far from being detectable, because of the technical difficulties and 
limitations of the time. Contemporary improvements in technology have made possible to propose new ideas to reveal 
the Lense-Thirring precession by analyzing the data-sets on the orbits of Earth satellites (see e.g. [l3| where, for the 
first time, the use of LAGEOS satellite was proposed). Several proposals have been recently published to measure 
the Lense-Thirring effect by natural and artificial bodies in some Solar System scenarios. For example, in the 
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Sun with Venus is considered. Mars with MGS spacecraft is discussed in [IH, while Jupiter with the Galilean moons 
(which is the original idea by Lense and Thirring) is studied in [TB|- Regarding the Earth with the existing LAGEOS 
and LAGEOSII satellites, recent results are reported in [13], while for the approved experiment LARES the expected 
forthcoming measurement are discussed in [la ]. 

On the other hand, the experiment Gravity Probe-B [llj has been devoted to another gravitomagnetic effect due 
to Earth's rotation, i.e. the Pugh-Schiff effect consisting of the precessions of the spins of four gyroscopes carried 
onboard t he sp caecraft [13, HH- This experiment has detected the effect and its magnitude in the gravitational field of 
the Earth [22|. The originally expected accuracy was 1% or better, but it is still unclear if it will be finally obtained 
because of unexpected systematic effects arisen in the data analysis. Other experiments (like GP-C [23l. [24 [25]) have 
been proposed to reveal the space-time structure, which is affected by gravitomagnetism, for example evidencing clock 
effects around a spinning massive object. In particular, concerning the so-called gravitomagnetic clock effect, we have 
to stress that its most investigated form consists of the difference between the orbital periods of two counter-rotating 
satellites. 

Recently, gravitomagnetic effects have been considered also in the framework of gravitational lensing. By using 
the Fermat principle and the standard theory of gravitational lensing, the gravitomagnetic corrections to the time 
delay function and the deflection angle for a geometrically thin lens can be derived. Such corrections can induce 
observational effects both in point-like [2(| and in extended gravitational lenses (as the isothermal sphere and the disk 
of spiral galaxies [13, Hll ■ Other researches concerning the gravitomagnetic effects on time delay and light deflection 
have been pursued. In [2§], the gravitomagnetic effects in the propagation of electromagnetic waves in variable 
gravitational fields of arbitrary- moving and spinning bodies have been studied, while, in [30l. l3ll [32I] . the gravitational 
lensing due to stars with angular momentum, and then inducing gravitomagnetic effects, have been considered. 

Finally, the analogy between general relativity and electromagnetism suggests that there is also a galvano- 
gravitomagnetic effect, which is the gravitational analog of the Hall effect. This effect takes place when a current 
carrying conductor is placed in a gravitomagnetic field and the conduction electrons moving inside the conductor are 
deflected transversally with respect to the current flow. Such a galvano-gravitomagnetic effect, considering current 
carrying conductors, could be used for detecting the gravitomagnetic field of the Earth. A discussion of the effect and 
its measurability is in [33l . [34], [35| . 

In this paper, we want to study how the relativistic theory of orbits for massive point-like objects is affected by 
gravitomagnetic corrections. In other words, we want to consider the orbital effects of higher-order terms in v/c and 
this is the main difference with respect to the standard gravitomagnetic effect so far considered. In this case, the 
problem of gravitomagnetic vector potential entering into the off-diagonal components gm of the metric g M „ can be 
greatly simplified and the corrections can be seen as further powers in the expansion in c _1 (up to c -3 ). Nevertheless, 
the effects on the orbit behavior are interesting and involve not only the precession at peri-astron but also nutation 
corrections as we will show below. This means that it could be misleading to neglect such effects when the weak 
field approximation is not so weak, as in the case of point-like compact objects moving in a tight-binding regime 
or spiralizing each other as in the case of evolved binary systems constituted by black holes and/or neutron stars. 
A study in this sense is in [36] where the possibility of measuring the Lense-Thirring effect with the double pulsar 
J0737-3039A is discussed. 

In particular, we can study the evolution of compact binary systems in the extreme mass ratio limit, i.e. the 
mass of the moving particle is m and the mass that produces the gravitational field is M, so that m < M. This 
constraint is satisfied by several real systems. For example, there has been gathering evidence suggesting the existence 
of supermassive black hole with masses in the ranges 10 6 -j- 1O 9 M ) in galactic nuclei [37ll3a]. One expects that small 
compact objects (1 -j- 10M©) from the surrounding stellar population will be captured by these black holes following 
many-body scattering interactions at a relatively high rate [33, S3] ■ 

Our approach suggests that, in the weak field approximation, when considering higher order corrections in the 
motion equations, the gravitomagnetic effects can be particularly significant, also in a rough approximation, giving 
rise also to chaotic behaviors in the transient regime dividing stable from unstable orbits. Generally, such contributions 
are discarded since they are assumed too small but they have to be taken into account as soon as the v/c is not so 
small. 

Sec. II is devoted to the discussion of the gravitomagnetic corrections which have to be considered when relevant 
mass-energy current effects are presented in a given problem. The geodesies, and then their spatial components, the 
trajectories, are corrected by such terms. We derive the Christoffel symbols with gravitomagnetic corrections and the 
vector form of geodesies. In particular, the metric "gravitomagnetically" corrected is achieved and the conditions in 
which the vector potential V 1 can be substituted with its point-like counterpart &v l /c where $ is the static Newton 
potential and v l the velocity of the test-particle m moving around the generator of the gravitational field M. 

In Sec. Ill, the theory of orbits is discussed. We review the Newtonian and the relativistic theory considering, in 
particular, the role of relativistic corrections [13, S3]- In Sec. IV, after constructing an effective Lagrangian coming 
from the line element with the gravitomagnetic effect, we derive the equations of motion. Numerical results for orbits 
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and their phase-space portrait are presented in Sec.V. Discussion and conclusions are drawn in Sec. VI. 

II. GRAVITOMAGNETIC EFFECTS 

Before treating the theory of the orbits with the gravitomagnetic effects, let us get some insight into gravitomag- 
netism and show how to derive the corrected metric. A recent book concerning both theoretical and experimental 
aspects of gravitomagnetism is (HI, while the Lense-Thirring effect is discussed in 

A remark is in order at this point: any theory combining, in a consistent way, Newtonian gravity together with 
Lorentz invariance has to include a gravitomagnetic field generated by the mass-energy currents. This is the case, of 
course, of General Relativity: it was shown by Lense and Thirring [lfl, Sfl, S3, Sa| , that a rotating mass generates 
a gravitomagnetic field, which, in turn, causes a precession of planetary orbits. In the framework of the linearized 
weak-field and slow-motion approximation of General Relativity, the ensemble of the so-called gravitomagnetic effects 
are induced by the off-diagonal components of the space-time metric tensor which are proportional to the components 
of the matter-energy current density of the source. It is possible to take into account two types of mass-energy currents 
in gravity. The former is induced by the matter source rotation around its center of mass: it generates the intrinsic 
gravitomagnetic field which is closely related to the angular momentum (spin) of the rotating body. The latter is due 
to the translational motion of the source: it is responsible of the extrinsic gravitomagnetic field. This concept has 
been discussed in Refs.fHol. IHlj. Then, starting from the Einstein field equations in the weak field approximation one 
obtain the gravitoelectromagnetic equations and then the corrections in the metric. Let us start from the weak field 
approximation of the gravitational field 1 

9nv(x) = + VW' \h^v{x)\ « 1. (1) 

where r/^ is the Minkowski metric tensor and << 1 is a small deviation from it [52]. 

The stress-energy tensor for perfect - fluid matter is given by 

T" u = (p + pc 2 ) u^u v - pg» u (2) 
which, in the weak field approximation p <C pc 2 , is 

T 00 ~ pc 2 , T 0j ~pcv j , T"~pi;V. (3) 
From the Einstein field equations — (8ttG / c 4 )T^ , one finds 

V ' h 00 = — r p, (4) 



8ttG 



V hij = —2-Sijp, (5) 



2 . IQitG , 
V h o 3 = — 5 flP v . (6) 

where V 2 is the standard Laplacian operator defined on the flat spacetime. To achieve Eqs. (@])-([6]), the harmonic 
condition 



f v K» = o » (7) 



Notation: latin indices run from 1 to 3, while greek indices run from to 3; the flat spacetime metric tensor is r/ M „ = diag(l, — 1, — 1, — 1). 
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has been used. 

By integrating Eqs. (H])-©, one obtains 



2$ 

"00 — 2 ' (°) 



2$ 

ftij = --^-Sij , (9) 



h 0j = ^5 jt V l . (10) 



The metric is determined by the gravitational Newtonian potential 



&(x) = —G I 9 n d 3 x' , (11) 



" = ~ G I ■ (12 » 

given by the matter current density pv l of the moving bodies. This last potential gives rise to the gravitomagnetic 
corrections. 

From Eqs(U) and l(8|l- fT2|) . the metric tensor in terms of Newton and gravitomagnetic potentials is 

ds 2 = U + |f) c2dt2 - ^^cdtdxi - (l - ^ 5i j dx i dx j . (13) 

From Eq.lfl3|) it is possible to construct a variational principle from which the geodesic equation follows. Then we 
can derive the orbital equations. As standard, we have 

x a + r^afx" = , (14) 

where the dot indicates the differentiation with respect to the afRne parameter. In order to put in evidence the grav- 
itomagnetic contributions, let us explicitly calculate the Christoffel symbols at lower orders. By some straightforward 
calculations, one gets 

1 oo — u 

r°- 



and by the vector potential V , 
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, dV j ^ 


c 3 V dxi 


"T" dx 1 J 


1 _9*_ 




dx k 




2 / dV k _ 


_ dV j \ 


c 3 \ dxi 


dx k J 


If 9* ! 


;k , 9* 


c* \dx~T L 


i "f" ~dx~i 



(15) 

r fe If 9* rfe I ./$ rfe _ 9* r \ 

y ~ ^ \dx~T u i f dx^ u j dx^ u V) 



In the approximation which we are going to consider, we are retaining terms up to the orders <£>/c 2 and /c 3 . 
It is important to point out that we are discarding terms like (§/c 4 )d§/dx k , (V J /c 5 )d$/dx k , (&/c 5 )dV k /dx J , 
(V k j c e )dV : > I dx 1 and of higher orders. Our aim is to show that, in several cases like in tight binary stars, it is not 
correct to discard higher order terms in v/c since physically interesting effects could come out. 
The geodesic equations up to c -3 corrections are then 

? d 2 t 2 9$ dt dx j 2 dV m „ dV rn \ dx z dx j 



da 2 c 2 dxi da da c 3 \ dxi dx 1 I da da 



~, 5 ™^- + S i™^-r -J--J- = , (16) 
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for the time component, and 



d 2 x k 1 » / dt\ 2 19$. dx i dx j 



da 2 c 2 dxi y da J c 2 dx k 13 da da 

2 d<i> dx l dx k t 4 ( dV k dV m \ dt dx i _ 



c 2 dx l da da c 3 \ dx J dx k J da da 

for the spatial components. 

In the case of a null-geodesic, it is ds 2 = da 2 — 0. Eq. lfT3|) gives, up to the order c 

W l . , ( , 2$ 



-3 



C(it = -^r-dx +1 t" ) dleuclid , (18) 



where dl 2 uclid — Sijdx l dx^ is the Euclidean length interval. Squaring Eq. (|18)1 and keeping terms up to order c 3 , 
one finds 



C ' dt ' Z = ( 1 - 3f ) dl eucUd + —dx l dl eucUd . (19) 



Inserting Eq. lfT9|) into Eq. (fT7l) . one gets, for the spatial components, 

d 2 x k 2 <9$ fdl euclld \ 2 2d<f>dx l dx k 4 f dV k c 8V m \ dl eucHd dx j 



6 m — = Q . (20) 

der 2 ' c 2 eta* \ da J c 2 dx l da da c 3 V cto Jm <9x fe J da da 

Such an equation can be seen as a differential equation for dx k jda which is the tangent 3- vector to the trajectory. 
On the other hand, Eq. (|20|) can be expressed in terms of leuciid considered as a parameter. In fact, for null geodesies 

and taking into account the lowest order in v/c, da is proportional to dl eU ciid- From Eq. i|16p multiplied for ^1 H — > 

we have 

1$^ _ ± 6im v m — 
da \ da c 2 da c 4 lm da 

and then 



-\- + -^~-^r n V m — )=0, (21) 



where, as standard, we have defined the affine parameter so that the integration constant is equal to 1 Substituting 
Eq. (|18)1 into Eq. (|22|) . at lowest order in v/c, we find 

dleuclid -. /OQ^ 

~cla~ ~ ' {l6) 
In the weak field regime, the spatial 3-vector, tangent to a given trajectory, can be expressed as 

dx k cdx k , , 

(24) 



da dleudid 
By defining 



Eq.j2Qj becomes 



de k ^d£_2_d£ ;fe ±fdXl_c dV r - ^ _ . 



dleuclid c 2 <9a; fe c 2 9x z c 3 V eta-? 9x fc 
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which can be expressed in a vector form as 

de 2 4 
=- [ V $-e(e-V$)] + -3 [eA(VAV)] . (27) 

The gravitomagnetic term is the second one in Eq . (|27)1 and it is usually discarded since considered not relevant. This 
is not true if v/c is quite large as in the cases of tight binary systems or point masses approaching to black holes. 
Our task is now to achieve explicitly the trajectories, in particular the orbits, corrected by such effects. 



III. THEORY OF ORBITS 

Orbits with gravitomagnetic effects can be obtained starting from the classical Newtonian theory and then correcting 
it by successive relativistic terms. Here we give, for the sake of completeness, a quick review of classical and relativistic 
theory of orbits showing how gravitomagnetic effects are the further corrections to be taken into account. A detailed 
discussion of classical and relativistic theory of orbits can be found in jH, 

A. The Newtonian theory 

The motion of a test particle in a spherically symmetric Newtonian gravitational field, can be achieved starting 
from a variational principle where the Lagrangian is 

„ 1 2 GM , . 

c =r+— (28) 

where the particle mass has been assumed unitary. The velocity, in spherical coordinates, is 

v 2 = r 2 + r 2 9 2 + r 2 sin 2 6ip 2 . (29) 

Here the dot denotes the ordinary derivatives with respect to the time. The Euler-Lagrange equations are easily 
derived. For 0- component, we have 

— (r 2 f?) = r 2 sin cos 0^> 2 , (30) 

where an obvious solution is = n/2; in fact the motion is plane and the variable cannot be taken in consideration 
any more. The equation 

! (•»#)-<>. (31) 

gives 

r 2 if = const = H , (32) 
which is nothing else but the conservation of the angular momentum. Finally, we have 

r = rip 2 =- . (33) 



It is convenient to introduce the new variable 



u{<p) = ~ (34) 
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Being 



and using Eq. ([34|) and Eq. ([32|) . it results 



, du , 



1 du 2 du dip 



-r 



2™"^ _ 2 • / 



u dt dtp dt 

From this equation, one gets 



r 2 ipu' = -Hu'. (36) 



d ( du\ TT d<p d ( du\ . ,, H 2 

r = —n 

and then Eq.((33]) is 



u +u = — (38) 



where the trivial solution u — (r = oo) is discarded. The solution of Eq. (|38j) is 



GM , 

u= -jjY + Bcos(tp - tpo), (39) 



and then, imposing tpo = 0, one gets the orbits in polar coordinates 



"foO = i ■ (40) 



GM 



Here k = ^ 2 and e is the ellipticity whose value can give elliptic, hyperbolic and parabolic orbits (55j. Summarizing 

the solution for 9 gives the planar motion, the solution for <p gives the angular momentum conservation, while the 
solution for r gives the orbits. 

B. The relativistic theory 

The relativistic case can be seen as a correction to the Newtonian theory of orbits. As before, we can start from a 
Lagrangian which can be deduced from the Schwarzschild line element, that is 

C = e v (i ) 2 - e A (rf - r 2 (() 2 + sin 2 . (41) 
The Euler-Lagrange equation for 9 is 

4- (r 2 9) = r 2 sin 6 cos 9ip 2 . (42) 
as V / 

In analogy with Eq. (|30|) (the two equations differ for ds in place of dt), the solution of this equation is 9 = 7r/2; 
again, as in the classical case, the motion is plane and 9 disappears as dynamical variable. The equations for x° = ct 
and x 3 = <p admit two first integrals of motion since the Lagrangian does not depend on x° and on x 3 but only on 
their derivatives. We have 

R \ 

l- — \±o=l, r 2 ip = h, (43) 
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corresponding to the first integrals of energy and angular momentum. R s is the Schwarzschild radius. For 
x 1 = r we can use the definition C = g^ v x^x v = 1 instead of the corresponding second order equation. Being 

e v = e~ x = ( 1 ) , we have 



r 



C=[l-^) (x ) 2 - . W ' -r 2 (0 2 +sin 2 fly> 2 ) =1. (44) 



Replacing Eq. ((43j) and considering = 7r/2, we have 

«■-<•-$(>-*)-(-*)• 

As in the Newtonian case, using the variable given by Eq. ([34|) and using the second of Eqs. (|43| . it is 

f = —hu . (46) 



Inserting Eq.(@6]) and Eq.((34]) in Eq.(|i5]). we get 



V - bfv! - h z u z (1 - R s u) = (1 - R s u) . (47) 

This equation gives, by a quadrature, the solution u = u(ip) with the periastron precession but, in order to compare 
the result with the Newtonian case, we can derive Eq. (|47|) considering that f = —hu 2 u". One obtains 

u" + u=^- 2 +\r s u* (48) 
This equation can be easily compared with the corresponding Newtonian case ([38)1 since 

1 H 

h~r 2 -tp=—. (49) 

c c 

2GM 

Being Rs = — 5 — , it follows that 



R s _ (<^£\ (<? \ _ (GM\ 
2h 2 ~ \ c? ) \H 2 J ~ { c 2 ) 



(50) 



This means that the relativistic correction to the test particle motion is due to the second member of (|4"8)l . Such a 
term is small is small if compared to the other. In fact, using (|49|l we have 



2 R SU 2 _ _. 2 2 3H 



2 



J|| r 2 c 2 



= 3h 2 W~^^ = 3 - , (51) 



so we can use a perturbation approach to deal with it. As said, such a relativistic correction is responsible for the 
perihelion precession. However, in strong field and high relative velocity regime, such term has relevant effects. 



C. Relativistic corrections due to gravitomagnetic effects 

Starting from the above considerations, we can see how gravitomagnetic corrections affect the problem or orbits. 
Essentially, they act as a further v/c correction leading to take into account terms up to c -3 , as shown in Sec. II. 
Let us start from the line element (fT3| which can be written in spherical coordinates. Here we assume the motion 

GM 

of point-like bodies and then we can work in the simplified hypothesis $ = and V 1 = &v l . It is 
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</»•'-' = ( 1 + ^p-^ cdt 2 - ^1 - ^p-^j [dr 2 + r 2 d9 2 + r 2 sin 2 9 dtp 2 ] - ^ cdt {[cos 6 + sin 9 (cos (p + sin (p)] dr 
[cos (cos y> + sin (/?) — sin 9] rd9 + [sin 9 (cos y> — sin p)) rdp] 



As in the Newtonian and relativistic cases, from the line element ([52|) . we can construct the Lagrangian 



, 2$\ / 2$ 

1+ -5- t 2 - ( 1 



c 2 J \ c 2 



f 2 + r 2 9 2 + r~ 



sin 2 0(^ 2 — { [cos + sin 9 (cos 95 + sin ip)\ \ 



+ [cos 9 (cos + sin p) — sin 0] r9 + [sin (cos <p — sin 99)] r<p j . (52) 

. . / 2$\ 

Using the relations (|43[) and being, as above, C = 1, one can multiply both members for I 1 + — J . In the planar 

motion condition 9 — ir/2 , we obtain 

r - [ 1 + — J I 1 2") h + -j J - — [(cos V + snip) r - (cos p-smp)p}= I 1 + J , (53) 

and then, being — ^- = and u = - it IS 

c r r 

; 2 - /i 2 (l - i? 2 u 2 ) (V 2 + w 2 ) + 4R ^ Ul [(cos 93 + sin p) u + (cos p> - sin u 2 ] = (1 - R s u) . (54) 

By deriving such an equation, it is easy to show that, if the relativistic and gravitomagnetic terms are discarded, the 
Newtonian theory is recovered, being 

u" + u = ^. (55) 

This result probes the self-consistency of the problem. However, it is nothing else but a particular case since we have 
assumed the planar motion. This planarity condition does not hold in general if gravitomagnetic corrections are taken 
into account. 

IV. ORBITS WITH GRAVITOMAGNETIC EFFECTS 

From the above Lagrangian (|52| , it is straightforward to derive the equations of motion 



c (rc 2 + 



GM) (J) 2 + sin 2 6></> 2 ) r 2 (56) 



c cos Or 



cr (rc 2 + 2GM) 

-AG Mi ((cos6>(cos0 + sin</>) - sin 9)9 + sin9(cos(f> - sin r + cGM f 2 - cGAli 2 



(c cot 9 (rc 2 + 2GM) 9<t>r 2 + r ( 2GM esc #(sin — cos <p)i + cr (rc 2 + GM) <\> 

r 2 (rc 3 + 2GMc) ' 

(rc 2 + 2GM) sin 9(j) 2 + r UgM(cos 6»(cos + sin 0) - sin 6)i - 2cr (rc 2 + GM) &\ 



r 2 (rc 3 + 2GMc) 



(57) 



(58) 



corresponding to the spatial components of the geodesic Eq. 1(20]) . Due to the numerical calculations which we are 
going to perform below, we consider the explicit form of the equations of motion. We have not considered the time 
component t since it is not necessary for the discussion of orbital motion. 
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As remarked above, from C = 1 the first integral f is achieved. It is: 

1 



± 



rl64G 4 M 4 r 



r 2 c 6 + AG 2 M 2 (-c 2 + 4 sin 20(cos + sin 0) + 4 sin 2 sin 20 + 4) 
((2 cos 20(cos + sin 0) + sin 26* sin 20)0 + (2 cos 20 sin 2 + sin 20(cos - sin 0)) 
(r 2 c 6 + 4G 2 M 2 ( - c 2 + 4sin20(cos0 + sin0) + 4 sin 2 sin 20 + 4)) (r 3 (0 2 + sin 2 00 2 )c 6 - 4GMc 4 + 
-r((E 2 -2)c 6 + 4G 2 Af 2 ((c 2 + 4 sin 20 (cos + sin 0) - 4 cos 2 sin 20 - 4^0 2 - 8sin0(cos0 - sin0) 

)S 0(cos + sin - sin 0)00 + sin 2 0(c 2 + 4 sin 20 - 4) 2 )) ^ - 8G 2 M 2 r ((2 cos 20(cos + sin 0) + 

1 
2 

sin 20 sin 20)0 + (2cos20sin 2 + sin20(cos0 - sin0))0^) , (59) 

which is the natural constrain equation related to the energy. The double sign comes out from the quadratic form of 
the Lagrangian. For our purpose, the positive sign can be retained. 

In the following calculations, we adopt geometrized units. Our aim is to study how gravitomagnetic effects modify 
the orbital shapes and what are the parameters determining the stability of the problem. As we will see, the energy 
and the mass, essentially, determine the stability. Beside the standard periastron precession of General Relativity, a 
nutation effect is genuinely induced by gravitomagnetism and stability greatly depends on it. A fundamental issue 
for this study is to achieve the orbital phase space portrait. 

V. NUMERICAL RESULTS 

The solution of the above system of differential equations presents some difficulties since the equations are stiff and 
their numerical solutions can diverge in several test points. Some numerical algorithms allow to change dynamically 
the meshing in order to decrease the mesh size near the critical points. 

For our purposes, we have found solutions by using the so called Stiffness Switching Method to provide an automatic 
tool of switching between a non-stiff and a stiff solver coupled with a more conventional explicit Runge-Kutta method 
for the non-stiff part of our differential equations. 

We have used for the computation the 6 th version of Wolfram Software Mathematica package jHsJ. The stiffness of 
the differential equations is evident from Fig. [TJ where the first and second derivative of r, plotted with respect to t, 
show steep peaks corresponding to the points where the radial velocity changes its sign abruptly. We show the time 
series of both r(t) and r(t) together with the phase portrait r = f(r) and r(t), assuming given initial values for the 
angular precession and nutation velocities (see also Fig. [5l In FigHJ the results for a given value of nutation angular 
velocity with a time span of 10000 steps is shown. It is interesting to see that, by increasing the initial nutation 
angular velocity, being fixed all the other initial conditions, we get curves with decreasing frequencies for r(t) and 
r(t). This fact is relevant to have an insight on the orbital motion stability (see Fig|4|). We have taken into account 
the effect of gravitomagnetic terms, in Fig. [2j showing the basic orbits (left) and the orbit with the associated velocity 
field in false colors (right). From a rapid inspection of the right panel, it is clear the sudden changes of velocity 
direction induced by the gravitomagnetic effects. 

To show the orbital velocity field, we have performed a rotation and a projection of the orbits along the axes of 
maximal energy. In other words, by a Singular Value Decomposition of the de-trended positions and velocities, we 
have selected only the eigenvectors corresponding to the largest eigenvalues and, of course, those representing the 
highest energy components (see Fig [21). 

The above differential equations for the parametric orbital motion are non- linear and with time- varying coefficients. 
In order to have a well-posed Cauchy problem, we have to define: 

• the initial and final boundary condition problems; 

• the stability and the dynamical equilibrium of solutions. 

We can start by solving the Cauchy problem, as in the classical case, for the initial condition putting f = , = 0, 
= and = f and the result we get is that the orbit is not planar being 0^0. In this case, we are compelled to 
solve numerically the system of second order differential equations and to treat carefully the initial conditions, taking 
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into account the high non-linearity of the system. A similar discussion, but for different problems, can be found in 

[H3,[5i|. 

A series of numerical trials on the orbital parameters can be done in order to get an empirical insight on the orbit 
stability. The parameters involved in this analysis are the mass, the energy, the orbital radius, the initial values of 
r, 4>, 9 and the angular precession and nutation velocities <j> and 9 respectively. We have empirically assumed initial 
conditions on r, <f> and 9. 

The trials we have performed can be organized in two series, i.e. constant mass and energy variation and constant 
energy and mass variation. 

• In the first class of trials, we assume the mass equal to M — IMq and the energy E n (in mass units) varying 
step by step. The initial orbital radius tq can be changed, according to the step in energy: this allow to find out 
numerically the dynamical equilibrium of the orbit. We have also chosen, as varying parameters, the ratios of 
the precession angular velocity </> to the radial angular velocity r and the ratio of the nutation angular velocity 
9 and the precession angular velocity 4>. The initial condition on 4> nas been assumed to be = and the 
initial condition on 9 has been Qq = | . For M = 1 (in Solar masses) , t = 5 and <j> = — jq, we have found out 

two different empirical linear equations, according to the different values of 9, <j>. We obtain a rough guess of 
the initial distance r = r (E n ) around which is possible to find a guess on the equilibrium of the initial radius, 
followed by trials and errors procedure. 

• In the second class of trials, we have assumed the variation of the initial orbital radius for different values of 

f 

mass at a constant energy value equal to E n = 0.95 in mass units. With this conditions, we assume cf> — — 

and assume that 9 takes the two values 1/2 and 1/10. We can approach the problem also considering the 
mass parameterization, at a given fixed energy, to have an insight of the effect of mass variation on the initial 
conditions. The masses have been varied between 0.5 and 20 Solar masses and the distances have been found 
to vary according to the two 3rd-order polynomial functions, according to the different values of 9 with respect 
to the mass. 

In summary, the numerical calculations, if optimized, allow to put in evidence the specific contributions of gravit- 
omagnetic corrections on orbital motion. In particular, specific contributions due to nutation and precession emerge 
when higher order terms in v/c are considered. 

VI. DISCUSSIONS AND CONCLUSIONS 

In this paper, we have discussed the theory of orbits considering gravitomagnetic effects in the geodesic motion. In 
particular, we have considered the orbital effects of higher-order terms in v/c which is the main difference with respect 
to the standard approach to the gravitomagnetism. Such terms are often discarded but, as we have shown, they could 
give rise to interesting phenomena in tight binding systems as binary systems of evolved objects (neutron stars or 
black holes). They could be important for objects falling toward extremely massive black holes as those seated in the 
galactic centers [13, The leading parameter for such correction is the ratio v/c which, in several physical cases 
cannot be simply discarded. For a detailed discussion see for example [26l. I27I I2H l30| . A part the standard periastron 
precession effects, such terms induce nutations and are capable of affecting the stability basin of the orbital phase 
space. As shown, the global structure of such a basin is extremely sensitive to the initial angular velocities, the initial 
energy and mass conditions which can determine possible transitions to chaotic behaviors. Detailed studies on the 
transition to chaos could greatly aid in gravitational wave detections in order to determine the shape, the spectrum 
and the intensity of the waves (for a discussion see [HqI. IgoI] ) . 

In a forthcoming paper, we will discuss how gravitomagnetic effects could affect also the gravitational wave pro- 
duction in extreme gravitational field regimes. 
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Figure 1: Plots of r(t) and r(t) for a test mass M — IMq, energy per mass unit E n = 0.95 and initial values for the orbital 
radius ro = 20, given in terms of Schwarzschild radius. The initial values of the angular precession velocity cj> and the angular 
nutation velocity 8 have been chosen according to the following criterium: assuming a given value of the initial radial velocity 
r, the initial values of the angular precession velocity and of the angular nutation velocity are <j> = — j^f and 6 — —-^f = jQ<j)- 
The phase portrait of r = f(r) is shown. The adopted time span is 10000 steps. 
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Figure 3: Breaking points examples: on the left panel, the first four orbits in the phase plane are shown: the red one is labelled 
I, the green is II, the black is III and the fourth is IV. As it is possible to see, the orbits in the phase plane are not closed 
and they do not overlap at the orbital closure points; we have called this feature breaking points. In this dynamical situation, 
a small perturbation can lead the system to a transition to the chaos as prescribed by the Kolmogorov-Arnold-Moser (KAM) 
theorem [3]. On the right panel, it is shown the initial orbit with the initial (square) and final (circles) points marked in black. 
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Figure 4: Plots of orbits with various energy values. For each value of energy, four plots are shown: the first on the left column 
is the orbit, with the orbital velocity field in false colors. The color scale goes from blue to red in increasing velocity. The 
second on the left column is the orbit with a different nutation angular velocity. On the right column, the phase portraits 
r = r(r(t)) are shown. Energy varies from 0.3 to 0.4, in mass units. The stability of the system is highly sensitive either to very 
small variation of ro or to variation on the initial conditions on both precession and nutation angular velocities: it is sufficient 
a variation of few percent on ro to induce system instability 
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Figure 5: Time series r = r(t) and phase portrait f = /(r) (top panels), time series r — r(t) and r = r(t) (middle panels), with 
the 3D orbits (bottom panels) for a Chandrasekhar mass M = 1.4 in solar units. We assumed the following initial conditions: 
r'o = —1/10, 4>o = — fo/10 while we have performed two trials assuming, for the initial condition on the nutation angular velocity 
#o, two limit values which we have found, according to our empirical procedure, i.e. 9q = 0o/2O and 9o = <f>o/2 respectively. At 
the bottom, the 3D orbits are plotted (left panel with 8o = 4>o/l0 and the right panel with 8o = </>o/2.) 



